//#include <dynamic_matrix.hpp>
#include <matrix.hpp>

#include <iostream>
#include <iomanip>
#include <algorithm>
#include <iterator>
#include <boost/lambda/lambda.hpp>

#include <random.hpp>

using namespace std;
using namespace sm;
using namespace boost::lambda;

int main()
{
    cout << setprecision(5);

    dynamic_matrix<long double> a(11,11);

    cout << a << endl;
	cout << "\nsizeof dynamic_matrix<long double>(11,11) is " << sizeof(a) << endl;
    long double b[]= { 1, 2, 3, 4, 5, 6, 7, 8, 9 };

    copy( b, b+3, a.col_begin(1));

    copy( b+3, b+8, a.row_begin(4));
    
    dynamic_matrix<long double> c = a;
    //fill( c.begin(), c.end(), 100.0 );
    fill( c.rbegin(), c.rend(), 100.0 );
    copy( b+3, b+8, c.col_begin(4));

    dynamic_matrix<long double> d = -c;
    
    //cout << a << endl << c  << endl << d << endl;

    copy( d.upper_diag_begin(4), d.upper_diag_end(4), ostream_iterator<long double>(cout, "\t"));
    cout << endl;
    copy( d.lower_diag_begin(0), d.lower_diag_end(0), ostream_iterator<long double>(cout, "\t"));
    cout << endl;
   
    c[5][5]= 13.3;

    cout << c << endl;
    cout << endl;

    dynamic_matrix<int> e = c;
    cout << e << endl;
    cout << endl;
    
    d+=c;
    cout << d << endl;
    cout << endl;
   
    dynamic_matrix<long double> f = d + 10.0;

    cout << f << endl;

    dynamic_matrix<long double> g = f*d;

    cout << g << endl;
    cout << endl;

    dynamic_matrix<long double> h = g/87.0;

    h += 1;

    cout << h << endl;
    cout << endl;
    
    h.transpose();

    h -=1;
    
    cout << h << endl;


    dynamic_matrix<long double> i = h/g;

    fill( i.diag_begin(), i.diag_end(), 1.3 );

    dynamic_matrix<long double> j = 1.0 / i;

    cout << i << endl;
    cout << endl;
    
    cout << j << endl;
    cout << endl;

    dynamic_matrix<long double> k(17,17);
    fill( k.begin(), k.end(), 0 );

    dynamic_matrix<long double>l(k);
    fill( k.diag_begin(), k.diag_end(), 10 );
    fill( l.diag_begin(), l.diag_end(), 0.1 );

    cout << l/k << endl;
  
    k += 100;
    l += 100;

    fill( k.lower_diag_begin(2), k.lower_diag_end(2), 987.1203 );
    fill( k.col_begin(5), k.col_end(5), 0.0039848 ); 
    fill( k.row_begin(5), k.row_end(5), 9990.0039848 ); 

    cout << "k^(-1)=\n" << k.inverse() << endl;
    cout << "det(k)=\n" << k.det() << endl;
    cout << "l^(-1)=\n" << l.inverse() << endl;
    cout << "det(l)=\n" << l.det() << endl;
  
    dynamic_matrix<long double> o(k/l);

    cout << "k=\n" << k << endl;
    cout << "l=\n" << l << endl;
    cout << "o=k/l=\n " << k/l << endl;
    cout << "o*l = \n" << o*l << endl;

    dynamic_matrix<long double> m(17,17);
    fill( m.begin(), m.end(), 990 );
    cout << "m=\n" << m << endl;
    cout << "m^(-1)=\n" << m.inverse() << endl;
    cout << "det(m)=\n" << m.det() << endl;

    dynamic_matrix<long double> n(2, 2);
    fill( n.begin(), n.end(), 0 );
    fill( n.diag_begin(), n.diag_end(), 100 );
    n[0][0]= -1;
    cout << "n = \n" << n << endl;
    cout << "det(n) = \n" << n.det() << endl;

    dynamic_matrix<long double> p( 11, 11 );

	feng::variate_generator<long double> vg;
    fill( p.begin(), p.end(), vg() );
    fill( p.diag_begin(), p.diag_end(), 190 );
    fill( p.upper_diag_begin(7), p.upper_diag_end(7), 19 );
    fill( p.lower_diag_begin(4), p.lower_diag_end(4), 1909 );
    cout << "p=\n" << p << endl;
    cout << "p^(-1) = \n" << p.inverse() << endl;
    cout << "det(p) = \n" << p.det() << endl;

    dynamic_matrix<long double> p_inv = p.inverse();

    cout << setprecision(5);
    cout << "P P^(-1) = \n" << p * p_inv << endl;
	cout << "\np=\n" << p << endl;

	//dynamic_matrix<long double, 0> q(20, 20);
	//generate( q.begin(), q.end(), ran_gen() );
	//dynamic_matrix<long double> r = q.inverse();
	//cout << "q=\n" << q << endl;
	//cout << "q^-1=\n" << r << endl;
	//cout << "q*q^{-1} = \n" << q * r << endl;

	dynamic_matrix<int> s( 11, 11 );
	int index = 0;
	for_each( s.begin(), s.end(), _1 = var(index)++ );
	cout << "\ns=\n" << s << endl;

	dynamic_matrix<int> r(12,12);
	cout << "\nbefore copy, r=\n" << r;
	copy( s.row_range_begin(1,5)+3, s.row_range_end(1,5), r.row_range_begin(0,4)+1 );

	cout << "\ns=\n" << s << endl;
	cout << "\nrow [1,5]+3, {1,5} -> row [0,4]+1:: r=\n" << r << endl;

	dynamic_matrix<int> t(13,13);
	cout << "\nbefore copy, t=\n" << t;
	copy( s.col_range_begin(1,9)+5, s.col_range_end(1,9), t.col_range_begin(0,8)+2 );

	cout << "\ns=\n" << s << endl;
	cout << "\ncol [1,9]+5, {1,9}-> col [0,8]+2:: t=\n" << t << endl;

	dynamic_matrix<int> u(14,14);
	cout << "\nbefore copy, u=\n" << u;
	copy( s.row_range_begin(2,6)+5, s.row_range_end(2,6), u.col_range_begin(0,4)+3 );

	cout << "\ns=\n" << s << endl;
	cout << "\nrow [2,6]+5, {2,6} -> col [0,4]+3 :: u=\n" << u << endl;

	dynamic_matrix<int> v(15,15);
	cout << "\nbefore copy, v=\n" << v;
	copy( s.col_range_begin(4,10)+5, s.col_range_end(4,10), v.row_range_begin(0,6)+4 );

	cout << "\ns=\n" << s << endl;
	cout << "\ncol [4,10]+5, {4,10}-> row [0,6]+4:: v=\n" << v << endl;

#if 0
	cout << "\nreverse range iterator test \n";

	dynamic_matrix<int> w(15,15);
	cout << "\nbefore copy, w=\n" << w;
	copy( s.col_range_begin(3,9)+5, s.col_range_end(3,9), w.row_range_rbegin(1,7)+4 );
	cout << "\ns=\n" << s << endl;
	cout << "\ncol [3,9]+5, {3,9}-> row r[1,7]+4:: w=\n" << w << endl;

	dynamic_matrix<int> x(15,15);
	cout << "\nbefore copy, x=\n" << x;
	copy( s.col_range_rbegin(2,8)+5, s.col_range_rend(2,8), x.row_range_rbegin(0,6)+4 );
	cout << "\ns=\n" << s << endl;
	cout << "\ncol r[2,8]+5, r{2,8}-> row r[0,6]+4:: x=\n" << x << endl;

	dynamic_matrix<int> y(15,15);
	cout << "\nbefore copy, y=\n" << y;
	copy( s.col_range_begin(4,10)+5, s.col_range_end(4,10), y.col_range_rbegin(0,6)+4 );
	cout << "\ns=\n" << s << endl;
	cout << "\ncol  [4,10]+5,  {4,10}-> col r[0,6]+4:: y=\n" << y << endl;

	dynamic_matrix<int> z(15,15);
	cout << "\nbefore copy, z=\n" << z;
	copy( s.row_range_begin(4,10)+5, s.row_range_end(4,10), z.col_range_rbegin(0,6)+4 );
	cout << "\ns=\n" << s << endl;
	cout << "\nrow [4,10]+5, {4,10}-> col r[0,6]+4:: z=\n" << z << endl;

	dynamic_matrix<int> A(15,13);
	cout << "\nbefore copy, A=\n" << A;
	copy( s.row_range_rbegin(6,10)+5, s.row_range_rend(6,10), A.col_range_rbegin(2,6)+4 );
	cout << "\ns=\n" << s << endl;
	cout << "\nrow r [6,10]+5, {6,10}-> col r[2,6]+4:: A=\n" << A << endl;

#endif




    return 0;
}





